Ab initio theory of coherent phonon generation by laser excitation 
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We show that time-dependent density functional theory (TDDFT) is applicable to coherent optical 
phonon generation by intense laser pulses in solids. The two mechanisms invoked in phenomenolog- 
ical theories, namely impulsively stimulated Raman scattering and displacive excitation, are present 
in the TDDFT. Taking the example of crystalline Si, we find that the theory reproduces the phe- 
nomena observed experimentally: dependence on polarization, strong growth at the direct band 
gap, and the change of phase from below to above the band gap. We conclude that the TDDFT 
offers a predictive ab initio framework to treat coherent optical phonon generation. 
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There has been much experimental progress in the 
study of intense electromagnetic fields interacting with 
condensed matter using pump-probe techniques on fem- 
tosecond time scales These interactions are a chal- 
lenging subject for theory, in view of the need to go be- 
yond pcrturbativc methods in dealing with strong fields. 
One promising theoretical approach useful to describe 
electron dynamics on femtosecond time scales is time- 
dependent density functional theory (TDDFT) Q. In 
this Letter we apply the TDDFT to the generation of co- 
herent phonons by strong laser pulses. Our goals are both 
to test the utility of the TDDFT in this domain and to 
assess the validity of phenomenological models that are 
in current use. In the past, two mechanisms have been in- 
voked to explain the generation of coherent phonon 0, 0] . 
The impulsively stimulated Raman scattering was pro- 
posed for the coherent phonon generation in dielectrics 
with a laser pulse whose frequency is lower than the di- 
rect band gap. In this mechanism, electrons are virtu- 
ally excited following adiabatically the laser electric field. 
The crucial quantity is the Raman tensor, the derivative 
of dielectric function with respect to the phonon coordi- 
nate. The other mechanism, called displacive excitation, 
requires higher frequencies to generate real electron-hole 
excitations in the final state [5|-|7[ . These excitations then 
shift the equilibrium position of the phonon coordinates. 
In this work we consider a bulk Si under irradiation of 
laser pulses of frequencies below and above the direct 
band gap, and show that the TDDFT is computationally 
feasible, includes two above-mentioned mechanisms, and 
produces results that are in qualitative agreement with 



experiments [H The TDDFT calculations prove to 
be also useful to evaluate phenomenological and macro- 
scopic models for the phonon generation process. 

Our computational framework is based on equations 
of motion derived from a Lagrangian for a periodic crys- 
talline system under a time-dependent, spatially uniform 
electric field [lOj . The Lagrangian is 



L = 



df{(en ion - en e )4> - E xc [n e ]} 



dR a 
dt 



r * — * 



dR a 7* 

e — — A . 

dt 



(1) 



Here ipi is the time-dependent electron orbitals, taken 
as Bloch orbitals in a unit cell of volume J7. n e (r,t) — 
J2i \i>i(r,t)\ 2 represents the electron density distribution. 
R a are atomic positions. The electromagnetic field terms 
are split into a long-range spatially uniform part A(t) 
and a periodic part given by a Coulomb potential </>. 
Variations with respect to the orbitals ipi, potential </>, 
and atomic coordinates R a result in the time-dependent 
Kohn-Sham equation for ipi, the Poisson equation for </>, 
and the Newton equation for R a , respectively. All the 
equations except those for R a are the same as those em- 
ployed in and Q. 
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FIG. 1: Geometry of the electric field and the optical phonon 
displacement in the 8-atom unit cell. The [Oil] x [100] plane 
and atoms on the plane are drawn with small arrows which 
show the direction of the optical phonon coordinate. 



To introduce the external laser field, we express the 
vector potential A(t) as a sum of an external field A ext (£) 
and the induced field A in< i(t), with A{t) = A cxt (t) + 
And(i) and treat A- m ^(i) as dynamic. The variation with 
respect to A- m d(t) yields the following equation of motion, 



» d 2 A ind (t) 
4ttc 2 dt 2 



f dr\j lon -j e }~^N e A(t) (2) 



e 

c ./,, 



To simulate the time-dependent electric field of the laser 
pulse, we take A oxt (t) to have the form 



A. 



dt'Sg sin 2 (~^r ) sin^'/' 



(3) 



for < t < T p and zero otherwise, 



with T p = 16 fs and 



So corresponding to peak intensity I = 10 12 W/cm 2 . 

The laser pulse is directed on a [100] Si surface at nor- 
mal incidence with a linear polarization oriented along 
the [Oil] axis. We show in Fig. [T]the atomic positions 
of Si atoms in the plane defined by the [Oil] and [100] 
axes. The 4 atoms lying on the plane are shown. The 
optical phonon coordinate which couples to the laser field 
is shown by vertical blue arrows. 

Our calculations are based on the LDA density func- 
tional 12[, treating the four valence electrons of Silicon 
explicitly and using the Troullier-Martins pseudopoten- 
tial 13[ . We employ the real-time and real-space scheme 
which was developed by us [l4| . The geometry is taken to 
be a simple cubic unit cell containing 8 Si atoms, with lat- 
tice constant a = 10.26 au. We have carefully examined 
the convergence of the results with respect to numerical 
parameters. We find that a spatial division of 16 3 , k- 
space grid of 24 3 , and the time step of At = 0.08 au is 
adequate for our purposes, and these numerical parame- 
ters are adopted for the results reported below. To make 
the present calculation feasible, parallel computation dis- 
tributing /c-points into processors is indispensable. We 
note the calculated direct band gap of Si is 2.4 eV, smaller 
than the measured value of 3.3 eV. 
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FIG. 2: Electric fields are shown as a function of time. The 
red solid line shows the applied laser pulse, Eq. ([3]), char- 
acterized by the peak intensity, I = 10 12 W/cm 2 , frequency 
huj — 2.5 eV, and the pulse duration, T v = 16 fs. The green 
dashed line shows the summed electric field of applied and 
induced ones, multiplied by a factor 15. 
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FIG. 3: Top panel shows the ground-state electron density 
in the plane shown in Fig. [T] The middle and bottom pan- 
els show the change of the electron density from that in the 
ground state by the laser pulse described in Fig. [2] The mid- 
dle panel corresponds to the time A and the bottom panel to 
the time B in Fig. [5] respectively. In the middle and bot- 
tom panels, the red color indicates the increase of the electron 
density, while blue color indicates the decrease. 
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We first show the electron dynamics induced by a laser 
pulse. Figure [2] shows the time dependence of the elec- 
tric fields. The red solid curve shows the electric field 
of applied laser pulse E e7Ck (t) = —(l/c)dA QXt /dt. We 
choose the laser frequency hu — 2.5 eV, close to the 
value of the direct band gap. The green dashed curve 
shows the sum of the applied and induced electric fields, 
Etot{t) = E ext (t) + E in d(t). The difference of the magni- 
tudes of the two fields comes from a dielectric screening. 

Figure |3] shows the electron density in the plane of Fig. 
[TJ The top panel shows the ground-state electron den- 
sity, and the middle and bottom panels show the change 
of electron density from that in the ground state at two 
times, marked A and B in Fig. [2j respectively. In the 
middle and bottom panels, red and blue indicate an in- 
crease or decrease of electron density, respectively. At 
time A, the electric field is maximum and there is a strong 
virtual excitation of the electrons. In the middle panel 
of Fig. [3l a movement of electrons is seen in the bond 
connecting two Si atoms. At the time B, the external 
electric field ended. Since the ultrashort laser pulse in- 
cludes frequency components above the direct band gap, 
there appear real electron-hole excitations. In the bottom 
panel of Fig. [3l one can see that the excitation results in 
a decreased density in the bond region and an increase 
near the Si atoms but away from the bond. One should 
note that the coloring of the middle and bottom figures 
are different by a factor of 40 to improve the visibility of 
the density change at time B. 

We next examine how the character of the electronic 
excitation changes as the laser frequency increases from 
below to above the direct band gap. Characteristics of 
the excitation as a function of time are shown for frequen- 
cies huj = 2.25 eV, 2.5 eV, and 2.75 eV in Fig. gj The 
top panel shows the total increase in energy in the unit 
cell, including both electronic excitation energy and the 
electromagnetic field energy. The red solid curve shows 
the results for a frequency below the band gap. Here the 
energy drops almost to zero after the pulse is over, as to 
be expected. The green dashed curve, corresponding to 
a frequency at the band gap, shows that some excitation 
energy remains after the end of the pulse, comparable in 
magnitude to the total energy at the peak. Finally, the 
blue dotted curve shows that above the gap the laser- 
electron interaction is highly dissipative, leaving a large 
excitation energy in the final state. The middle panel 
in the figure shows the number of excited electrons as a 
function of time. This is calculated by taking the over- 
laps of the time-dependent occupied orbitals with the 
initial state static orbitals as in Ref. 11] . The results 



are qualitatively very similar to what we found for the 
energy. Below the direct band gap, the excited electron 
shows a peak during the pulse and then drops off to a 
very small value in the final state. At higher frequencies, 
the excitations remain in the final state and it is not pos- 
sible to distinguish the real excitation from the virtual 



0.016 
0.012 
0.008 
0.004 



ons 


o 


o 


0.006 












"O 






0.004 


u 




X 





o 


0.002 



JD 




E 









z 


0.02 


nm] 


0.015 


> 






0.01 







o 




o 


0.005 


LL 








(a) 



2.25 eV 
2.5 eV 
2.75 eV 



1 H^H 1 (- 



(b) 



(c) 



5 10 15 20 25 30 
Time [fs] 

FIG. 4: Electron excitation of the crystal during and after the 
pulse for several laser frequencies across the direct band gap. 
The top panel (a) shows the energy in the unit cell including 
electron-hole excitation energy and the electric field energy. 
The middle panel (b) shows the the number of electron-hole 
pairs in the unit cell. The bottom panel (c) shows the force 
on the optical phonon coordinate. 



one during the pulse. In summary, one sees an adiabatic 
response below the gap switching rather abruptly to a 
strongly dissipative response above the gap. 

Finally, in the bottom figure, we show the calculated 
induced force for the three frequencies. Note that the 
ion positions are fixed in these calculations; the accelera- 
tions are small and the resulting displacements would be 
inconsequential. The lowest frequency, shown by the red 
solid curve, gives a force envelope that follows the shape 
of the pulse intensity. This is just what one would expect 
from the adiabatic formula [3j, Eq. (2)]. One also sees 
high frequency oscillations superimposed on the envelope 
of the curve. The frequency of these oscillations are twice 
the laser frequency, again as expected from the adiabatic 
formula. The green dashed curve shows the force for a 
laser frequency of fiuj — 2.5, nearly at the direct band 
gap. One still sees a large peak at 10 fs associated with 
instantaneous high field intensity. However, there is a 
residual force after the end of the pulse which is rather 
constant with time. This is just what one expects for dis- 
placive mechanism. At this point, we have shown that 
TDDFT reproduces at a qualitative level the role of the 
two mechanisms. Beyond that, the relative sign associ- 
ated with them can be extracted from the graph. The 
last case shown, hui = 2.75, is 0.35 eV above the direct 
gap. Here the displacive mechanism is completely dom- 
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FIG. 5: The amplitude (a) and the phase (b) of the phonon 
oscillation Eq. as a function of laser frequency cj. 

inant, although one can still see an enhancement of the 
force during the pulse. 

We now integrate the time-dependent force to get the 
lattice distortion associated with the phonon coordinate. 
In principle, the restoring potential for the lattice vibra- 
tion is included in the evolution equations, but the am- 
plitude of the lattice displacement is too small numeri- 
cally to include it in the direct integration. So for this 
part of the analysis we simply assume a harmonic restor- 
ing potential consistent with the observed optical phonon 
frequency, f ph onon = 15.3 THz. 

To analyze the characteristics of the coherent phonon, 
we fit the oscillation of the displacement to a cosine func- 
tion as in conventional parametrization of the experimen- 
tal reflectivity measurements 0, 

q(t) = -qo cos(w p/l i + 0) +q. (4) 

Fig. [5] shows the amplitude and phase as a function 
of laser frequency, fitted in the time interval 40-90 fs. 
Below the direct gap energy the phase is close to it/ 2 
as expected for the Raman mechanism. The amplitude 
remains almost constant in this frequency region, also 
consistent with the Raman mechanism. One sees a quite 
sharp drop from that value to </> = as the direct gap is 
crossed, showing the transition to the displacive behav- 
ior. The amplitude also shows a sudden increase across 
the direct gap. Several experimental measurements are 
also shown on the figure for the phase. Two of them 
[l,[l5[ are in the Raman regime. The theory supports the 
results of Ref . Q , which reports a value close ton/2. The 



other measurement does not appear consistent with our 
theory or indeed with the other experiment. The phase 
has also been measured in the gap region 0, shown by 
the square on Fig. [5jb). This point should be compared 
with the theory at the corresponding calculated gap en- 
ergy, 2.4 eV. In both theory and experiment the phase 
has decreased from the Raman value, but decrease seems 
larger for the experimental measurement. Both results 
are in a range where the mechanism is changing rapidly. 
All in all, we find the agreement quite satisfactory on a 
qualitative level, particularly since the phase could have 
come out with an opposite sign {(f) « it). 

At higher frequencies, the theoretical phase goes to 
zero as expected for the displacive mechanism, but then 
it rises again beyond 4.5 eV, approaching tt at hcu — 
5 eV. There is a corresponding dip and growth in the 
amplitudes associated with a change in the sign of the 
displacive force. Different electron orbitals are excited at 
the high frequency, and apparently those orbitals have 
an opposite sign contribution to the displacive shift. 

We also examined the dependence of amplitude and 
phase of the coherent phonon on the intensity of the laser 
pulse. At all frequencies we examined, the amplitude of 
the phonon is proportional to the laser intensity. In the 
impulsive Raman mechanism which is applicable below 
the band gap, this dependence is expected from the adia- 
batic formula [3j, Eq. (2)]. In the displacive mechanism, it 
is also expected if the medium is excited by a one-photon 
absorption process. The phase of the coherent phonon is 
found to be sensitive only on the frequency but not on 
the intensity until the multiphoton absorption processes 
become significant. 

In summary, we have derived and carried out a compu- 
tational method to apply time-dependent density func- 
tional theory to laser-lattice interactions, taking as an 
example the excitation of coherent optical phonon by 
femtosecond-scale laser pulse in silicon. The qualitative 
agreement between theory and measured phase of the 
coherent phonon confirms the utility of the TDDFT to 
describe electron dynamics resulting from intense laser 
pulse in solids. 

The numerical calculation were performed on the 
massively parallel cluster T2K-Tsukuba, University of 
Tsukuba, and the supercomputer at the Institute of Solid 
State Physics, University of Tokyo. TO acknowledges 
support by the Grant-in-Aid for Scientific Research No. 
21740303. GFB acknowledges support by the National 
Science Foundation under Grant PHY-0835543 and by 
the DOE under grant DE-FG02-00ER41132. 
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